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ABSTRACT 

We model three dSph galaxies, Draco, Ursa Minor and Fornax, as axisym- 
metric stellar systems embedded in spherical dark-matter potentials, which are 
in dynamical equilibrium without significant external tidal forces. We construct 
non-parametric two- and three-integral models that match the observed surface- 
density profiles and the current kinematic samples of ~ 150—200 stars per galaxy; 
these models naturally produce the so-called "extra-tidal extensions" , which had 
previously been suggested as evidence of tidal stripping. Isochrone, NFW and 
power-law models fit the data, but we strongly rule out any centrally condensed 
mass distribution like a dominant central black hole, as well as constant-density 
and (for Draco and Fornax) constant mass-to-light ratio models. The average V- 
band mass-to-light ratio is 400 ± 80 Mq/Lv,© within 0.75 kpc for Draco, 5801}qo 
M Q /Ly )Q within 1.1 kpc for Ursa Minor and 25^ M Q /Ly )0 within 2.5 kpc for 
Fornax. Two-integral models fit the data almost as well as three-integral models; 
in contrast to previous suggestions we do not find that anisotropy contributes 
substantially to the high mass-to-light ratios in these dSph galaxies. 

Subject headings: Methods: statistical - galaxies: kinematics and dynamics 
- galaxies: fundamental parameters - galaxies: individual (Draco, Ursa Minor, 
Fornax) - cosmology: dark matter 



1. Introduction 

Dwarf galaxies include dwarf irregular (dlrr) and dwarf spheroidal (dSph) galaxies. 
Their luminosities (10 3 — 1O 8 L , see e.g. Belokurov et al. 2007) are much smaller than the 
characteristic luminosity L* ~ 2 x 10 10 L Q for giant galaxies. The faintest dSph galaxies 
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have luminosities comparable to globular clusters, but are more extended (half-light radius 
~ 0.5 kpc compared to ~ 0.005 kpc for globular clusters) and thus have much lower surface 
brightness. 

The formation history of dwarf galaxies is not well-understood. The standard theory 
is that dwarf galaxies formed at the centers of subhalos orbiting within the halos of giant 
galaxies. However, the standard ACDM model predicts a much larger number of subhalos 
than the observed number of dwarf galaxies. There are ~ 50 known dwarf galaxies in 
the Local Group, compared to 300 predicted subhalos with circular velocity ~ 10 — 20 
kms-^Klypin et al. 1999). 

Aaronson (1983) was the first to determine the velocity dispersion of a dSph galaxy, 
using three stars in the Draco dwarf (6.5 kms" 1 ). Aaronson's result was later confirmed by 
more extensive observations (Olszewski, Pryor & Armandroff 1996; Wilkinson et al. 2004). 
It is found that all dwarf galaxies have central velocity dispersions ~ 6 — 25 km s _1 (Mateo 
1998). These dispersions are believed to reflect motions of the stars in the gravitational 
potential of the galaxy rather than (say) atmospheric motions or orbital motion due to 
a binary companion (Olszewski, Pryor & Armandroff 1996; Mateo 1998). If the systems 
are in dynamical equilibrium, the mass derived from these velocity dispersions is much 
larger than the stellar mass, resulting in a V-band mass-to-light ratio ranging from 7 to 500 
M /L y>o (Irwin & Hatzidimitriou 1995; Wilkinson et al. 2002; Palma et al. 2003; Kleyna et 
al. 2005; Lokas et al. 2005). Therefore, dSph galaxies are probably the darkest objects ever 
observed in the universe, and thus provide unique probes of the dark-matter distribution on 
small scales. 

Two possible alternatives to dark matter in dSph galaxies have been proposed. The 
first is Modified Newtonian Dynamics (MOND, Milgrom 1983). However, MOND has diffi- 
culty reproducing the different mass-to-light ratios of different dSph galaxies, mostly because 
MOND is based on a characteristic acceleration whereas the dSph galaxies seem to have a 
characteristic velocity dispersion instead. Moreover, MOND is difficult to reconcile with the 
power spectrum of fluctuations in the CMB (Spergel et al. 2006). The second alternative is 
tidal forces from the host galaxies. Tides have changed the structure and internal kinematics 
of the Sagittarius dwarf galaxy (Ibata, Gilmore & Irwin 1994; Ibata et al. 2001). So it is 
natural to postulate that other dSph galaxies are also undergoing tidal disruption, which 
leads to the high velocity dispersions (Kroupa 1997; Martmez-Delgado et al. 2001; Gomez- 
Flechoso & Martmez-Delgado 2003). However, (i) Sagittarius is much closer to the center of 
the Milky Way than other dwarf galaxies; (ii) tides should introduce strong velocity shear 
and dSph galaxies rotate slowly if at all; (iii) no mechanism of increasing velocity dispersions 
by tides appears to work for all dwarf galaxies (Mateo 1998; Klessen et al. 2003). There are 
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a few pieces of putative evidence that tides may influence the structure of Ursa Minor and 
Draco ( Martinez- Delgado et al. 2001; Palma et al. 2003; Muhoz et al. 2005). For exam- 
ple, "extra-tidal extensions" (e.g., Irwin & Hatzidimitriou 1995; Wilkinson et al. 2004) are 
found beyond a "tidal radius" by fitting a King profile (King 1966) to stellar number density 
profiles. However, the King profile was intended to fit relaxed, constant mass-to-light ratio 
systems such as globular clusters (Trager et al. 1995), and perhaps the profile can not be 
applied to dSph galaxies. 

In this paper we model the mass distribution of three dSph galaxies, specifically, Draco, 
Ursa Minor and Fornax. Since we focus on the inner parts (< 2 kpc) which are very regular 
(Irwin & Hatzidimitriou 1995) and not significantly affected by tidal stripping (Read et 
al. 2006), we assume that the inner part is in dynamical equilibrium with no significant 
external tidal forces, and test whether this assumption is self-consistent with a given mass 
distribution. 

A wide variety of approaches have been used to model the kinematics of dSph galaxies. 
Different analysis techniques and assumptions have led to different results. Taking Draco 
as an example (see also Figure 4), if we calibrate the results assuming a V-band luminosity 
of 1.8 x 1O 5 L (Irwin & Hatzidimitriou 1995), Odenkirchen et al. (2001) derive a V-band 
mass-to-light ratio of Fy= 195 ± 56 M /Ly i0 within 1.2 kpc using a King model. Wilkinson 
et al. (2002) and Kleyna et al. (2002) assume parametric forms for both potential and 
distribution (hereafter DF), and model Draco as a near-isothermal sphere (v C i rc oc r 017 ), 
with fy = 440 ± 240 M /Ly i0 within 0.7 kpc, more than a factor of two larger than derived 
by Odenkirchen et al. (2001) at a larger radius. By fitting velocity moments, Lokas (2002) 
investigates the dark matter distribution with parametric DFs and the assumption that dSph 
galaxies are spherical and derives f v — 175 ± 35 M /Ly i0 within 0.4 kpc. Mashchenko, Sills 
& Couchman (2004) propose that the mass of Draco could range from 7 x 10 7 to 5 x 1O 9 M , 
assuming an NFW (Navarro, Frenk & White 1997) potential, which has T v = 540 ± 170 
and 510 ± 160 at r = 0.75 and 1.1 kpc, respectively. Lokas et al. (2005) parametrize 
anisotropy by a constant anisotropy parameter f3 to fit velocity moments and found the 
mass-to-light ratio was ~ 370 M /Ly i0 and almost constant with radius for Draco. Wang et 
al. (2005) have determined dark matter distributions, using a sophisticated non-parametric 
mass distribution, assuming spherical symmetry. However, they assume that the distribution 
function is isotropic in velocity space, an assumption that can introduce much larger errors 
than the ones that are controlled by the use of non-parametric mass distribution. We believe 
that a better method is to use a non-parametric distribution function and a parametric 
potential. 
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The apparent ellipticities 1 of the stellar distributions of dSph galaxies vary from 0.13 
(Leo II) to 0.56 (Ursa Minor) (Mateo 1998). However, the gravitational potential is dom- 
inated by the dark halo in these galaxies, and there is no reason to suppose that the dark 
halo shape is the same as the shape of the stellar distribution. Hence for simplicity, we 
assume that the galactic potential is spherical. In addition, we assume that the distribution 
of stars in the dSph galaxy is axisymmetric. Other methods designed to model axisymmetric 
galaxies (e.g., Cretton et al. 1999) are more general in that they allow both the potential 
and the distribution function to be axisymmetric; however, these are designed to work for 
a continuous, rather than discrete, observed distribution of stars. We describe and test our 
method in §2 and §3, respectively. Then we derive the mass density profiles for Draco, Ursa 
Minor and Fornax in §4. We discuss the results in §5 and conclude in §6. 



We introduce a Cartesian coordinate system, with origin at the center of the dwarf 
galaxy and z— axis along the line of sight. The x and y-axes are the apparent major and 
minor axes in the plane of the sky, respectively, which can be determined from isodensity 
contours of the surface number density profile (isopleths). 

Then the distribution function may be written as f(E,L,L z >), where E is energy, L 
is scalar angular momentum, the z'— axis is the intrinsic symmetry axis of the stellar dis- 
tribution and L z i is the angular momentum along the z'-axis. Therefore, the z'-axis lies 
in the y — z plane. For an observational data set of positions and line-of-sight velocities 
{xi,yi,v Z i}, we model f(E,L,L z >) and determine both the parameters of the potential and 
the inclination 9, which is the angle between the z'-axis and the z-axis. 

Wu & Tremaine (2006, hereafter Paper I) have used maximum likelihood analysis to 
deal with the analogous problem for spherically symmetric DFs in a spherical potential, 
where the DF is a function of two integrals, E and L. It is straightforward to generalize this 
method to the case of an axisymmetric DF f(E, L, L z r). We define a probability function, 



which gives the probability density to observe a star at (x, y) with a line-of-sight velocity v z . 
Note that g is also a function of 9 and the assumed gravitational potential <3>. We maximize 



l li the apparent semi-major and minor axes of isopleths are a and b, then the apparent ellipticity e a = 
1 — b/a. Intrinsic ellipticity is defined similarly, i.e., it is equal to e a while the galaxy is seen edge-on. 



2. 



Method 





- 5 - 



the log likelihood function 

LH($, f,9) = J2 Vu /, 0), (2) 

i 

assuming parametrized potentials and non-parametric DFs as follows. 

For the potential, we assume an analytical form $(r, X) and infer the parameters X. 
We investigate six possible forms. 



1. An isochrone model, 



$(r) = GM , X = {M,6}, (3) 



where M is the total mass. 
2. The Navarro, Frenk & White (1997) density-potential pair, 

p(r) = - 



2 ' 

r I r ' 
— 1 + — 



$(r) = -47^^ V r rn/ , (4) 



r 

where r n is the concentration radius and X = {p„,r„}. 
3. A power-law density-potential pair 

p( r ) = Po( — 

2-Q 



•M = M ^ r °. i f I («<3), (5) 



(2 — a) (3 — a) \r 
where r is an arbitrary radius, and then X = {p , a}. 

4. A constant-density model, which is a special case of the power-law model with a = 
and X = {po}- 

5. A constant mass-to-light ratio model (see §4.3). 

6. A model in which the potential is dominated by a central black hole of mass M^h, 

$(r) = -GM bh /r, X={M bh }. (6) 
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For comparison with earlier analysis by other authors, we shall also fit the data using 
a more conventional method: fitting to the King (1966) three-parameter family of isotropic 
stellar systems, by minimizing x 2 f° r the surface-density and velocity-dispersion profiles (see 
§4.2). 

There is little rotation in dSph galaxies. Therefore, we do not attempt to model the 
observed rotation: we simply assume / is an even function of L z >. Furthermore, the non- 
rotating model has even symmetry around the x and y-axes, i.e., g(x,y,v z ) = g(\x\, \y\, \v z \), 
so we can assume that x, y, v z are all positive and < 6 < tt/2. These simplifications increase 
the computational efficiency significantly. 

To construct the DF, we divide the {E, L, \L Z //L\} space into Ne x Nl x N z > bins, which 
are denoted by the triple index mnl, m = 1,...,Ne, n = 1,...,Nl, I = l,...,N z i. The bins 
are mutually exclusive and cover all of the allowed E — L — L z i / L space, as well as some 
unallowed regions. In our calculation, bin mnl is defined by 

E m -i < E < E m , (7) 
^Ll(E m ) <L 2 < ^LKEJ, (8) 

> cos7/, where 7 Z = (9) 



cos7i_i > 



i z' 



Here E m _ 1 and E m are the lower and upper limits of the energy in the bin, and L C (E) is the 
maximum angular momentum at energy E, corresponding to a circular orbit. To achieve 
uniform spatial resolution for different potential models, we choose the boundaries of energy 
bins {E m } as the potentials at a set of evenly separated radii. The limits in L z //L are 
chosen so that the resolution in 7, the angle between the z' — axis and the vector normal to 
the orbital plane (the orbital inclination) is uniform. Note that for spherical systems, the 
natural choice is uniform in cos 7, not 7. However, when modeling flattened systems, this 
choice yields too few bins near 7 = 0. The factor u determines the maximum value for 7 
and must be 1 for the bins to cover all allowed regions. However, it can be chosen to be < 1 
to model a galaxy with very flat (disk-like) stellar part with higher computational efficiency 

Assuming that the DF is constant in bin mnl and zero elsewhere, we can calculate a 
probability distribution in (x,y,v z ) space (cf. Paper I), 

g mn i{x,y,v z \X,e) oc / dv x dv y dz. (10) 

Jhm. mnl 

we normalize each g mn i(x, y, i> z |X, 6) so that f g mn i(x,y,v z \X.,9)dxdydv z = 1. 
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Assuming the non-negative weights in the bins are W = {w mn i}, we have 
^2 Wmnl = 1 > 



m,n,l 



g(x,y,v z ) = ^2w mn ig mn i(x,y,v z ), 



({xi,yi,v zi }\X,6) 



(11) 



where LH is the log of the likelihood. Notice that g(x,y,v z ) has been normalized since all 
gmni{x,y,v z ) are normalized. For a fixed inclination 9 and fixed potential specified by the 
parameters X, we maximize LH(X, W, 9) with respect to W, then vary X and 9 to search 
for a global maximum LH(X,W,9). More generally, the posterior probability of (X, 9) is 
given by 



corresponding to an isotropic distribution of orientations. Our rule for choosing Px is to 
assign a uniform probability distribution for dimensionless parameters and a uniform proba- 
bility distribution in log scale for dimensional parameters. We choose Pw = 1 in this paper; 
that is, in contrast to Paper I we do not use regularization to smooth the DF (in Paper I, 
we found that the errors in derived potential parameters for samples of a few hundred stars 
were dominated by statistical noise rather than different regularization levels). 

It is straightforward to account for observational errors in velocities and observational 
cuts on the sample volume in velocity and position on the sky (velocity cuts are normally 
applied to eliminate interlopers). The former is dealt with by convolving g mn i(x,y,v z ) with 
a Gaussian error profile. The latter requires us to normalize g mn i(x,y,v z ) within the cuts in 
velocity and position. 

If the surface density of the kinematic sample is proportional to the overall surface 
density of stars, we use Equation (11). If not, the distribution of positions contains no 
information, so in the likelihood function we use the conditional probability 




(12) 



where Pg, Px and Pw are prior probabilities. We use 



Pe = sin#, 



(13) 



h(v zi x i: y h X,9) = — 

J g{xi,yi,v z \X,9)d 




(14) 



In this paper, we use Equation (11) for the calculations in this paper. 
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An observed surface number density profile may be available for many more stars, over 
a larger survey area than the kinematic sample. We can easily add this as an additional 
constraint. Assume the apparent ellipticity of the stellar part is e a , and that its surface 
number density profile along the x— axis is H b s (R e ) with error a^(R e ), where R e is the 
distance along the x— axis. We maximize the quantity 

Q = LH ~\x\ (15) 

Vz(Rej) 2 



x 2 



Ri&Re) 



Rp 



\J cos 2 + sin 2 0/(1 — e a ) 2 



wd\ f g(Rcos(f),Rsm(f),v z ) 

E ( Re > 00 / — 2a> ■ 2-777; ^dv z d(f) 

J cos z + sm 0/(1 — e a y 

Ef g mn i(R cos (j),R sin <f),v z ) 
w mn i / — — , 2 ^dv z d(p, (17) 
J cos 2 + sin 0/ 1 — e a Y 

m,n,l 

where R is the projected distance from the origin to the ellipse defining the isopleth at an 
angle relative to the x— axis. In Equation 17, S(i? e ) is free to be scaled by an arbitrary 
constant factor to minimize x 2 . For some systems (e.g., Draco), the available measured 
surface number density profile is averaged along concentric circles rather than isopleths, we 
need to simply set e a = in Equation (17). 



3. Simulations 

To test the algorithm in §2, we construct a dSph galaxy in a power-law potential with 

r = 0.5 kpc, 9 = = 0.7854, 
X = {p , a} = {8.5xlO 7 M /kpc 3 , 1.2}, 

f(E, L, L z .) = (V(£-™0) 2 /2002 + q ^-(£-1650)7200^ e sL,L c{ E) fW\ ^ (lg) 

where s is an anisotropy parameter and E is in units of (kms -1 ) 2 . The DF is chosen so that 
the energy distribution is bimodal to challenge our algorithm. 

Figure la shows the spatial distribution of 200 stars in the simulated galaxy with 
s = 5 (tangential anisotropy), which has e a ~ 0.3. We use the algorithm in §2 with 
(N E , N L , N z i) = (15,5,15) and no additional surface density data. Accounting for the 
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prior probability of 9 (eq. 13), 9 is derived to be 0.90 ± 0.07 by fitting a Gaussian to 
the probability distribution in Figure lb, which is consistent with the input value at the 
1.6(7 level. Figure lc shows the probability distribution of potential parameters, specifi- 
cally, X = {(8.1 ± 1.2)xlO 7 M /kpc 3 , 1.17 ± 0.25}. The triangle and plus sign are the input 
model and the best-fit model, respectively, which agree within the la error. Therefore, our 
algorithm can recover both potential parameters and inclination. 

Contrastingly, with the assumption that the stellar density is spherical and the DF 
is isotropic, i.e., (Ne, Nl, N z >) = (15,1,1), we may get large systematic errors. We have 
explored two cases. For the above input model with s = 5 (tangential anisotropy), the 
assumption leads to the best-fit model marked by the square with error bars in panel (d) 
of Figure 1. For another input model with s = — 5 (radial anisotropy) and the same input 
potential p , a, the best-fit model marked by the cross is also inconsistent with the input 
model. We conclude that three-integral models that allow for both velocity anisotropy 
and flattening must be used to avoid possible systematic errors in the derived potential 
parameters. 

For a face-on galaxy (e a = 0), we find that the 2a contour for potential parameters 
extends over at least two orders of magnitude in p . Because the line-of-sight velocity dis- 
persion is decoupled from the mass distribution for a face-on disk-like galaxy, an apparently 
round galaxy may be either a low-density spherical galaxy or a face-on high-density flat 
galaxy. Therefore, we get a poor constraint on potential parameters of an apparent round 
dSph galaxy. In practice, we never see disk-like dSph galaxies, so when modeling a galaxy 
with e a = we will not make a serious error by assuming that the stellar distribution is 
spherical and using the analogous methods described in Paper I for spherical systems of test 
particles 2 . This is not an issue for the three galaxies analyzed in this paper, which are suf- 
ficiently flattened that they yield well-constrained results for the potential parameters and 
inclination. 

4. Application to Draco, Ursa Minor and Fornax 

4.1. Data 

Our algorithm is applied to kinematic data {x, y, v z } and number density profiles £ ;> S (-R e ) 
of three dSph galaxies taken from the published literature. Table 1 lists several parameters 



2 On the other hand, for a flat edge-on galaxy (e a = 1.0), the intrinsic ellipticity et is also 1.0, so we may 
overestimate the mass by up to 50% if assuming that the stellar distribution is spherical (/ = f(E, L)). 
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of each galaxy, including the assumed heliocentric distance (D), the luminosity in the V- 
band (L v ), apparent ellipticity (e a ), and the position angle of the major axis derived from 
the isopleths, as usual measured eastward from north (PA). In the kinematic sample, we 
include stars, which have heliocentric velocities in the range v k ± Av and are inside an 
ellipse with semi-major axis Rk and ellipticity e a , to exclude possible interlopers which could 
strongly bias the results (the same cuts are applied to our likelihood distribution, so this 
procedure introduces no systematic errors). Table 1 also lists the number of stars (N s ) and 
annuli (N a ) in the surface number density profile £ of , s (i? e ), which has a maximum radius 
R m . For the outer parts (R > R g ), where the number density is very low, we combine a 
few annuli to increase the signal-to-noise ratio. For Draco, we exclude two bins in which the 
number density is more than 3a away from those in nearby bins. 



4.2. King model 

Since King models are spherical, we convert the elliptical surface density distribution to 
a spherical one, with a surface number density profile 

E'(R') = Z obs (R e ), R' = VT^~aRe, (19) 

so that the number of stars enclosed in a circle with radius R' and an ellipse with semi-major 
axis R are equal. A King model is defined by three parameters, which may be chosen as the 
characteristic velocity dispersion of the DF a, the difference between the potential at the 
tidal radius and that at the center $(0) in units of a 2 , and a core radius r c (see Binney & 
Tremaine 1987, for details). Since the DF is parametric, we do not use the method in §2, 
but fit velocity dispersion profiles and number density profiles, as shown in Figure 2. We 
have binned the kinematic data into iVj, = 10 bins with approximately the same number 
of stars in each bin. Considering the likelihood and assigning a flat prior probability for 
{<fr(0)/<r 2 , In a}, we estimate the parameters and their errors, which are given in Table 2. 

We fit the observed line-of-sight velocity dispersion profiles while Irwin & Hatzidimitriou 
(1995), having a smaller kinematic data set, only fit the core radial velocity dispersion. 
Moreover, Irwin & Hatzidimitriou (1995) estimate an average velocity dispersion 7.5 ± 1.0 
kms _1 from 46 stars for Ursa Minor, which is much lower than 13.3±0.4 kms _1 from 163 stars 
used in our calculation. Nonetheless, considering the factor of 2.0 in the velocity dispersion 
for Ursa Minor, the derive mass-to-light ratios (see Table 2) similar to the values of Irwin & 
Hatzidimitriou (1995), T v = 245 ± 155 M /Ly )O for Draco, 95 ±43 M Q /L V:Q for Ursa Minor 
and 7 ± 3 M /Ly i0 for Fornax. 

Figure 2 shows that a single-component King model does not fit the data of Draco 
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and Fornax, but is formally consistent with those of Ursa Minor. However, the last point 
of the Ursa Minor kinematic data suggests that a King model may be ruled out if future 
observations show that the flat dispersion profile extends to larger radii. In all three systems, 
the observational data have an excess in both velocity dispersion and surface number density 
at large radii, and the dispersion profile is flat rather than decreasing as in the King models. 



4.3. More general models 

We fit the parametrized models for the mass distribution in §2 to the data. We must 
assume a maximum radius R max , so that the upper energy limit of the bins in Equation (7) 
is $(Rmax)- We set R max as 2 kpc for Draco and Ursa Minor, 6 kpc for Fornax, which are 
reasonable because the surface number density is negligible at these radii. For the power-law 
model, we set tq as 0.28 kpc for Draco, 0.6 kpc for Ursa Minor and 1.1 kpc for Fornax to 
minimize the covariance of po an d ol in their bivariate probability distribution. 

The models labeled "constant mass-to-light ratio" require some explanation. Since our 
method assumes that the potential is spherical in this paper, we do not attempt to construct 
a full model of a non-spherical galaxy with constant mass-to-light ratio. To derive the 
potential of the constant mass-to-light-ratio model, we assume that all stars have the same 
mass-to-light ratio and construct a spherical density distribution as in Equation (19). So the 
luminosity profile v s {f) is proportional to the volume number density of stars. The spherical 
assumption is then used to derive the potential, but dropped when modeling the galaxy 
itself, i.e., we assume a non-spherical galaxy in a spherical potential. 

It is convenient to fit luminosity profile v s with a linear combination of two exponential 
functions 

( n P ~ br + C p~ dr < d 
Us{r) = \ r > R (20) 

where a, b, c, d are constant parameters listed in Table 3 so that 

L v = 4tt / u s (r)r 2 dr. (21) 
Jo 

The chi-square of the fitting shows that the luminosity profiles are indeed well approximated 
by Equation (20). 

This density distribution gives a potential 

• (r) = -i.GT v («Vr -<-*(* + 2M) + <Wr- t -y + 2M)\ 
\ tr a 6 J 
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The potential is dependent on X = {T v }, the constant mass-to-light ratio in the V-band. 

We consider two types of models, those with N L — 1, and N E , N z > 3> 1, which corre- 
spond to models having / = f(E,L z >) (two-integral models), and those with all of N L) N E) 
N z i 3> 1 (three-integral models). The two-integral models are much faster to compute and 
appear to fit the data just as well as the three-integral models for the three dSphs, so we 
focus on the two-integral models. 



4.4. Results 

4-4-1- Two integral models (Nl = 1) 

Table 4 lists the results for various models with (N E , N L , N z >) = (15, 1, 20), the difference 
in log likelihood ALH, potential parameters X, mass and mass-to-light ratio within r = 
(see Table 1). To determine the uncertainties in the derived parameters, we lay down points 
with a uniform probability distribution in ln(X) if X is dimensional and in X otherwise, 
then reject or save each point according to the likelihood at that point. From the saved 
points, we may estimate the model parameters and their errors. In practice, we consider a 
set of discrete values for 9, 

in 20 

9 = — , — arccos(e a ) < i < 10. (23) 

20 7T 

where i is a natural number. For all models of the three dwarf spheroidals, either i = 9 or 10 
gives the maximum P(X, 9) and others give much smaller likelihoods (3a away). Therefore, 
we do not try to fit 9 accurately, but give a rough estimate 9 ~ 1.45 ± 0.12. 

Several features of Table 4 are worth noting: (i) For the NFW and isochrone models, 
there is a strong correlation between the two parameters in X, which explains why the frac- 
tional error in M(r < Rk) is smaller than the fractional error in the potential parameters, (ii) 
The best-fit isochrone, NFW and power-law models have —4.2 < ALH < 0, corresponding 
to < 2.9o". The models of a dominant central black hole, constant density and (for Draco and 
Fornax) constant mass-to-light ratio are less likely at > 3.7a levels, (iii) The difference in 
log likelihood only indicates the relative likelihood among different models, so it is important 
to investigate whether the best-fit model is a good fit. In Figure 3 we show the observed 
velocity dispersion a Vt0 b s (R) in ten bins (top row). For the stars in a bin centered at R, 
say, {xj,yj,v Z j}, we calculate the predicted velocity dispersion a v j it (R) from the velocity 
distribution profile 

h(y z \R)=Y^g(x j ,y j ,v ie ). (24) 

3 
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Comparison of the first two rows of panels in Figures 2 and 3 shows that the present models 
provide a substantially better fit to the observed velocity-dispersion and surface- density 
profiles in all three galaxies than the best-fitting King models. This visual impression is 
confirmed in Table 6, which gives goodness-of-fit statistics for the velocity-dispersion profile 
(Xy) and surface-density profile (x 2 ) and their deviations from the expected values iVj, and 
N a (the number of bins), in units of the standard deviation \f2N~ b or y/2N a (la error). The 
deviations for the present models are < 2.1<r for all three galaxies, while the King models 
are much poorer fits for Draco and Fornax, (iv) Note that the models are fitted to the 
surface-density distribution as a function of radius, assuming the observed ellipticity e a (eqs. 
14-17). Thus it is not guaranteed that the the ellipticity of the isopleths in the models will 
be similar to the observed ellipticity (no doubt a better approach would have been to fit the 
two-dimensional surface- density distribution, but this is substantially more computationally 
expensive). Nevertheless the ellipticities of the models, quoted in Table 6 at the core radius, 
are similar to the observed ellipticities. 

The third row in Figure 3 shows the derived mass profiles of various models for the 
three dSph galaxies. All of the models shown are within 2. 5a of the best-fit models. The 
models have similar density distributions at intermediate radii, for example, 0.25 — 1.0 kpc for 
Draco, but (not surprisingly) differ substantially at radii beyond the outermost observational 
data. To summarize, our algorithm gives relatively small errors at intermediate radii, and 
the surface number density profile only weakly constrains the mass density profile beyond 
r = R k . 

The bottom row in Figure 3 shows the average mass-to-light ratio fV(r) for the best-fit 
models. The light inside a radius r is calculated from Equation (20). Comparing Table 2 
and 4, we can see that King model gives estimates of mass at the outer limit of the kinematic 
data, r = Rk, consistent with those from more general models within 2a. However, King 
model assumes a constant mass-to-light ratio, which is an average value. For the more 
general models, except at small radii where the mass density is very uncertain, the mass- 
to-light ratios increase with radius. Therefore, the global mass-to-light ratios are probably 
even greater than IV(-Rfc)- 

4 -4 -2. Three integral models (N L = 5) 

We have shown that two integral (Nl = 1) models can reproduce all of the observed 
features of dSph galaxies. We also list the results for three-integral models, (Ne, Nl, N z >) = 
(15,5,20), in Tables 5 and 6. The mass-to-light ratios are generally larger than the values 
for N L = 1 models by 10% to 50%. This result suggests that given the freedom to adjust 
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the DF in L direction, high mass models are better able to fit the observed data than low 
mass models, thus the derived mass distributions are more skewed to the high end. In the 
following sections, we adopt the conservative (lower) estimates of mass-to-light ratios that 
are derived with the case of N L = 1. 

5. Discussion 

We believe that the mass estimates presented here are the most accurate and reliable 
that have so far been obtained for the three dSph galaxies, Draco, Ursa Minor and Fornax. In 
our models, we assume that the stellar distributions are not significantly affected by the tidal 
force of the Milky Way. This assumption is confirmed by a dynamical argument: using the 
mass M(r < Rk) of the best-fit models in Table 4, which is a lower limit to the total mass, the 
Roche limit is 2.6 kpc (3.5Rt) for Draco, 2.6 kpc (2ARk) for Ursa Minor and 6.6 kpc (2ARk) 
for Fornax, assuming that that the Milky Way is an isothermal sphere with v circu i ar = 220 
kms -1 . The Roche limits also exceed R ma x, the maximum radius used in our dynamic 
models in §4.3. Therefore, our assumption that the dSph galaxies are isolated systems in 
virial equilibrium is self-consistent. Because our models can fit both radial velocity profiles 
and surface number density profiles, the so-called "extra-tidal extensions" in the surface 
number density profiles found by Irwin & Hatzidimitriou (1995) and Wilkinson et al. (2004) 
do not require any special ad hoc explanation. Thus it is not valid to consider them as the 
evidence of tidal stripping, as proposed by some authors (e.g., Martmez-Delgado et al. 2001; 
Gomez-Flechoso & Martmez-Delgado 2003; Muhoz et al. 2005). 

We strongly rule out any centrally condensed mass distribution like a dominant central 
black hole, as well as constant-density and (for Draco and Fornax) constant mass-to-light 
ratio models. For example, the constant-density model is ruled out at the level of 4.1a for 
Draco, 3.4<r for Ursa Minor and 4.4<r for Fornax. In contrast to the positive result of Lokas 
et al. (2005) based on the same data set, the constant mass-to-light ratio is ruled out at 8.6<r 
level for Draco. 

We derive an average mass-to-light ratio of 400 ±80 M /Ly iQ within 0.75 kpc for Draco, 
580l];oo M /Ly i0 within 1.1 kpc for Ursa Minor and 25^ M /Ly i0 within 2.5 kpc for Fornax. 
Note that these values are probably lower limits for the global mass-to-light ratio because 
the mass-to-light ratio increases with r for the best-fit models. The mass-to-light ratios we 
obtain for all three galaxies are larger than the estimates of global mass-to-light ratio by 
Irwin & Hatzidimitriou (1995) (245 ± 155 M /Ly i0 for Draco, 95 ± 43 M /Ly i0 for Ursa 
Minor and 7±3 M /Ly i0 for Fornax). The discrepancies arise because the results of Irwin & 
Hatzidimitriou (1995) and Odenkirchen et al. (2001) are based on a single-component King 
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model and because we have a larger sample of radial velocities, covering a larger range of 
radii. 

Figure 4 shows the derived mass-to-light ratio of Draco as a function of radius and the 
estimates by other authors. The mass-to-light ratio for Draco is assumed constant in some 
of these studies (Irwin & Hatzidimitriou 1995; Lokas et al. 2005) but is allowed to vary with 
radius in others (Odenkirchen et al. 2001; Lokas 2002; Kleyna et al. 2002; Mashchenko, Sills 
& Couchman 2004). Piatek et al. (2002) have proposed that the large mass-to-light ratios 
of dSph's may be partly due to velocity anisotropy, since anisotropy can inflate the derived 
Tyby a factor of three or so if the kinematic information is restricted to the central parts 
of the galaxy. The simulations of §3 confirm that velocity anisotropy and flattening of the 
stellar distribution (i.e., the dependence of the DF on L and L z >) must be accounted for 
to obtain accurate mass-to-light ratios. Nonetheless, for the three dSph galaxies, the DF 
f(E,L z t/L) fits the data reasonably well (§4.4), and the derived masses actually increase 
if we allow the dependence of DF on L (§4.3). So we find no evidence that three-integral 
models fit the current data significantly better than two-integral models, or that they prefer 
models with lower mass-to-light ratio. 

Wilkinson et al. (2004) and Kleyna et al. (2004) argue that the velocity dispersion 
profiles of Draco and Ursa Minor seem to show a sharp drop at radii of 0.81 and 0.84 kpc, 
respectively. The drop in the dispersion profile seems to coincide with a "break" in the 
number density profile, and Wilkinson et al. propose that this is evidence of a kinematically 
cold population. However, the sharp drop may only be present in particular binning schemes 
(Munoz et al. 2005) and is not seen with our binning scheme (Figure 2). Our models do not 
reproduce such a sharp drop. Thus, the drop is inconsistent with stationary axisymmetric 
models of the DF and also is not statistically significant given the present data. 

The large masses that we have derived may help alleviate but may not eliminate the 
problem that there are far fewer dwarf galaxies in the Local Group than predicted by standard 
cosmological models. 

6. Conclusion 

We model dSph galaxies as axisymmetric stellar systems in spherical potentials, which 
are in dynamical equilibrium without significant external tidal forces. For the three dSph 
galaxies, Draco, Ursa Minor and Fornax, plausible parametric models of the potential com- 
bined with non-parametric two-integral models of the stellar distribution function can nat- 
urally reproduce the observed velocity dispersion profiles and number density profiles. The 
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masses of these models are high enough so that tidal forces from the Milky Way do not sig- 
nificantly affect the observed kinematics. So the models are self-consistent and "extra-tidal 
extensions" simply reflect an extended axisymmetric envelope that is still in virial equilib- 
rium. The sharp drop in velocity dispersion which coincides with a "break" in the light 
distribution that has been claimed in the literature is not observed in our binning scheme 
and can not be reproduced in the scenario of virial equilibrium. 

The isochrone, NFW and power-law models all fit the data. We strongly rule out any 
centrally condensed mass distribution like a dominant central black hole, as well as constant- 
density and (for Draco and Fornax) constant mass-to-light ratio. The average V-band mass- 
to-light ratio in our best-fit two-integral models is 400 ± 80 Mq/Lv,© within 0.75 kpc for 
Draco, 580^}qo M /Ly i0 within 1.1 kpc for Ursa Minor and 25^ M /Ly iO within 2.5 kpc for 
Fornax. Our method considers both anisotropy of the DF and the non-spherical geometry 
of the stellar population, and investigates a large variety of mass distribution models, and 
thus should yield more reliable mass distributions than previous investigations. 

Our simulations show that anisotropy (dependence on L) and flattening (dependence on 
L z i) in the DF may substantially affect derived masses of axisymmetric galaxies. However, 
we detect no strong evidence of anisotropy in the three dSph galaxies, i.e., a two-integral DF 
f(E,L z i) is consistent with the data. Allowing anisotropy yields even larger mass-to-light 
ratios. 

Our method could be improved by fitting the two-dimensional surface number-density 
distribution rather than a one-dimensional profile. Observationally, it would be useful to 
extend the kinematic surveys to larger radii, since the current limits are still well within the 
radii where tidal forces are likely to complicate the analysis. 

I thank Scott Tremaine for numerous stimulating discussions and helpful comments on 
an earlier draft, James Gunn, David Spergel, Michael Strauss and Glenn van de ven for 
helpful discussions. I wish to thank Mark Wilkinson for providing the kinematic data of 
Draco and Ursa Minor and surface number density profile of Draco in tabular form. This 
research was supported in part by NASA grant NNG04GL47G and used computational 
facilities supported by NSF grant AST-0216105. 
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Fig. 1. — Algorithm tests, (a) A test dSph galaxy with 200 observed stars, described by 
Equation (18) with s = 5 (tangential anisotropy). (b) Probability distribution of 9 derived 
from the 200 stars. The input value is marked by the dotted line, (c) Likelihood contours 
for the potential parameters. The contours are na levels from the peak, i.e., e~ n / 2 of the 
maximum likelihood. The triangle is the input model and the plus sign is the best-fit model 
with (Ne, Nl, N z i) = (15,5,15). They agree within la error, (d) This panel shows the 
results of fitting flattened anisotropic galaxies with spherical, isotropic models. The triangle 
is the input model which is the same model described in (a) except that s can be +5 or 
—5 (radial anisotropy). The square with error bars is the best fit for the model with s — 5, 
assuming that the stellar part is spherical and the DF is isotropic, (N E , N L , N z >) = (15, 1, 1), 
while the cross is the same as the square except for the input model with s = —5. The large 
errors in the recovered potential parameters show that both anisotropy of the DF and non- 
spherical geometry of the stellar population should be accounted for whenever their effects 
may be significant. 
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Fig. 2. — Fitting King models to dSphs. The upper row shows the observed and best-fit 
velocity dispersion profiles a Vfi b s (points with error bars) and a v ju (solid lines). Also shown 
are the derived core radii (dashed lines) and those from Irwin & Hatzidimitriou (1995) (dash- 
dotted lines). The middle row shows the observed density profiles (E'(R'), eq. 19) and the 
best-fit. The bottom row shows the mass profiles of the best-fit King models. The core 
radius r c is 0.19 kpc for Draco, 0.15 kpc for Ursa Minor and 1.10 kpc for Fornax. A single- 
component King model does not fit the data of Draco and Fornax, but is consistent with 
those of Ursa Minor. Howeve, in all three galaxies the model shows a decline of velocity 
dispersion with radius that is not seen in the data. 
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R or r (kpc) 

Fig. 3. — Best-fit two integral models (N L = 1) for three dSph galaxies. Top row: observed 
velocity dispersion profiles a v ^ obs (R) and those of the best-fit models a v j it (R) (see eq. 24). 
Also shown are the core radius r c (dotted lines) and kinematic survey radius Rk (dashed 
lines). Second row: observed surface number density profiles £„;,<, (i? e ) and those of the 
best-fit models. Third row: derived mass profiles for the isochrone (solid lines), NFW 
(dotted lines) and power-law models (dashed lines). All three models agree relatively well 
at intermediate radii. Bottom row: derived average mass-to-light ratio within radius r. For 
Fornax, fyhas been multiplied by a constant 10.0 for clear presentation. The mass-to-light 
ratios tend to increase with r except at small radii, where the mass density is very uncertain. 
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Fig. 4. — Derived mass-to-light ratio for the best-fit isochrone (solid line), NFW (dotted 
line) and power-law (dashed line) models for two-integral models (Nl = 1) for Draco. Also 
plotted are measurements by a: Lokas (2002) b: Kleyna et al. (2002) c: Mashchenko, Sills & 
Couchman (2004) d: Lokas et al. (2005) e: Irwin & Hatzidimitriou (1995), f: Odenkirchen 
et al. (2001). The average mass-to-light ratio r is assumed constant over r in references d 
and e but is an estimate at the radii shown for other measurements. 



-23- 



Table 1. Data summary 11 





Draco 


Ursa Minor 


Fornax 


Reference 


galaxy parameters 


D (kpc) 


82 


66 


138 


1 


L v (1O 5 L V , ) 


1.8 


2.0 


140 


1,2 


e a 


0.29 


0.56 


0.30 


2 


PA (degree) 


88 


56 


41 


2,3 


kinematic data 


N k 


197 


163 


147 


4,5 


v k ± Av (kms -1 ) 


-290.7±30.0 


-245.2±35.0 


53.3±30.0 




Rk (kpc) 


0.75 


1.1 


2.5 




surface number density profiles 


N s 


18000 


1500 


42000 


2,4 


R m (kpc) 


1.2 


1.3 


4.0 




R g (kpc) 


n/a 


0.5 


2.5 




N a 


44 


38 


66 





a See section §4.1 for explanations of the parameters. 



References. - (1) Mateo (1998); (2) Irwin & Hatzidimitriou (1995); (3) Odenkirchen et 
al. (2001); (4) Wilkinson et al. (2004); (5) Walker et al. (2006). 



Table 2. Parameters of King models for dSph galaxies 



Galaxy 


*(0) 


a 


T . 


r t 


M tota i 


M/L 


M(r < R k 






(kms -1 ) 


(kpc) 


(kpc) 


(1O 6 M ) 


(M /L ) 


(1O 6 M ) 


Draco 


4.3 ±0.1 


11 ± 1 


0.182 ±0.005 


1.42 ±0.04 


34 ±5 


180 ± 30 


28 ±4 


Ursa Minor 


5.0 ±0.6 


18 ±2 


0.18 ±0.02 


1.7 ±0.3 


100 ±20 


500 ± 150 


83 ±28 


Fornax 


1.5 ±0.1 


23 ±2 


1.11 ±0.05 


2.84 ±0.04 


180 ±20 


11±4 


160 ± 50 
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Table 3. Parameters of luminosity profiles 





Rmax 


a 


b 


c 


d 


N a 


x 2 






(kpc) 


(L kpc- 3 ) 


(kpc- 1 ) 


(L kpc- 3 ) 


(kpc- 1 ) 






W2N a ) 


Draco 


2.0 


4.69 x 10 6 


9.17 


3.26 x 10 3 


0.98 


44 


68 


2.6 


Ursa Minor 


2.0 


6.27 x 10 6 


10.4 


1.01 x 10 4 


1.22 


38 


31 


0.8 


Fornax 3, 


6.0 


-4.21 x 10 7 


4.79 


4.21 x 10 7 


3.55 


66 


95 


2.5 



a We perform the optimization with the constraint that the central density a + c > 0. 
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Table 4. Results from two-integral models of dwarf spheroidal galaxies 



Mass model a 


ALH 


X{l} b 


X{2} c 


M(r < R k ) d 


T v (r < R k ) e 










(1O 6 M ) 


(M /L y>0 ) 


Draco 


1 


-0.3 


180+^° 


o.48i8;g 




4803?° 


2 


0.0 


1613 


2.ol? : g 


60 ± 12 


400 ± 80 


3 f 


-0.1 


77 ±8 


i-i^J 


79+24 
' z -16 


4803?° 


4 


-8.0 


80l£ 




77 +47 
' '-27 


510IS8 


5 


-13.7 


420 ± 60 




63 ±9 


420 ± 60 


6 


-36.8 


38 ±4 




38 ±4 


250 ± 30 


Ursa Minor 


1 


0.0 


14013 


U.18_ 004 


93^3 


pron+140 
OoU_ 1QO 


2 


-3.3 


9701™ 


20+ - 13 

u - zu -0.05 


93l?J 


cron+130 
°° u -100 


3 


-4.2 


22 ±4 


1 8+ - 2 
J-O-0.4 


1013? 


6303° 


4 


-20.0 


22+} 1 




lOOll 


630+ 340 

UOU_ 230 


5 


-4.0 


720 ± 60 




150 ± 10 


720 ± 60 


6 


-14.9 


80 ±9 




80 ±9 


500 ± 60 


Fornax 


1 


0.0 


660+ 370 


9 +0 - 2 

u - y -0.3 


340+?° 


253 


2 


-0.7 


13l^ 7 


o 0+1.4 
z - z -1.0 


3501™ 


253 


3 


-1.3 


19 ±2 


1 4+0- 2 
i - 4 -0.4 


3903o° 


283 


4 


-9.7 


20^ 




5001«8 


361S 


5 


-6.9 


17 ±3 




240 ± 40 


17 ±3 


6 


-17.7 


180 ± 30 




180 ± 30 


13 ±2 



a Model: 1. isochrone, 2. NFW, 3. power-law, 4. constant density, 5. constant mass-to- 
light ratio, 6. black hole. See §2 for description of the models. 

b The first parameter of mass distribution models: mass for isochrone and black hole models 
in units of 10 6 M , density for NFW, power-law, and constant density models in units of 10 6 
M kpc -3 , mass-to-light ratio for constant mass-to-light ratio model in units of M /Ly i0 . 
See Equations (3)- (6) for definitions of the parameters. 
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c The second parameter of mass distribution models: radius for isochrone and NFW models 
in units of kpc, density exponent for power-law model. 

d Mass interior to to the limiting radius R k of kinematic data. R k = 0.75 kpc for Draco, 
1.1 kpc for Ursa Minor and 2.5 kpc for Fornax. 

e The average mass-to- light ratio within r = R k . 

f The arbitrary parameter r (eq. 5) is chosen as 0.28 kpc for Draco, 0.6 kpc for Ursa 
Minor and 1.1 kpc for Fornax. 
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Table 5. Results from three-integral models of dwarf spheroidal galaxies 3. 



Mass model 


ALH 


X{1} 


X{2} 


M(r < i2 fc ) 


f v (r < R k ) 












(1O 6 M ) 


(M /L y , ) 


Draco 


1 


-0.3 


ooq+180 
oou -160 


0.51 


+0.08 
-0.16 




49011?° 


2 


-0.3 


2Qtf 7 


1.2; 


+2.0 
-0.5 


5 8 ;?° 


3801$° 


3 


0.0 


81 ± 14 


1.1 ±0.4 


77 ±26 


510 ± 170 


4 


-9.5 


oo+40 






oo+40 
° d -30 


rin+310 
Oiu -180 


5 


-5.2 


370 ± 40 






56 ±6 


370 ± 40 


6 


-31.9 


29 ±3 






29 ±3 


190 ± 20 


Ursa Minor 


1 


-0.3 


170 ± 50 


0.21 


+0.04 
-0.08 


105±» 


66011" 


2 


-0.1 


300™ 


0.37 


+0.12 
-0.19 




710118 


3 


0.0 


30^ 


1.4 ±0.3 


125l*> 


780118 


4 


-16.4 








100I8 


630+ 350 

UOU -230 


5 


-1.7 


720 ± 60 






115 ± 10 


720 ± 60 


6 


-11.1 


88 ±9 






88 ±9 


550 ± 60 


Fornax 


1 


-2.5 


2100+ 600 

ziuu -500 


1.5; 


f0.2 
-0.3 


5201^00 


38 ±7 


2 


-0.1 


c q+o.6 


4.0 ±0.4 


46018 


331^ 


3 


0.0 


25+? 


i q+0.2 
i -°-0.3 


5001™ 


36 + | 


4 


-13.0 


221J 1 






55011° 


401?^ 


5 


-9.7 


15 ±2 






210 ±30 


15 ±2 


6 


-17.8 


180 ± 40 






180 ±40 


13 ±3 



definitions and units of all parameters are the same as those in Table 4 
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Table 6. Statistics of the best-fit models 



Model 


e a 






\xl ~ N b \ 


N a 


x 2 


\X 2 ~ N a \ 










W2N b ) 






(V2Ay 


Draco 


King 


0.0 


10 


54.3 


9.9 


44 


108.5 


6.9 


f(E,L z ,) 


0.25 


10 


21.1 


2.5 


44 


38.7 


0.6 


f(E,L,L z .) 


0.30 


10 


21.0 


2.4 


44 


32.2 


1.3 


Ursa Minor 


King 


0.0 


10 


13.9 


0.9 


38 


40.3 


0.3 


f(E,L z ,) 


0.35 


10 


6.3 


0.8 


38 


53.2 


1.7 


f{E,L,L z ,) 


0.40 


10 


3.0 


1.6 


38 


52.1 


1.6 


Fornax 


King 


0.0 


10 


57.9 


10.7 


66 


172.5 


9.3 


f(E,L z .) 


0.15 


10 


13.1 


0.7 


66 


39.9 


2.3 


f(E,L,L z ,) 


0.30 


10 


12.9 


0.6 


66 


34.9 


2.7 



